MAE 271: Lab 2 - Inverted pendulum with vertical base excitation¶

Cooper Cook & Joshua Booth¶

Problem Introduction¶

The simulation problem below presents a vertical inverted pendulum with an excited base. This problem becomes interesting due to the fact that given certain initial angles from vertical and the right amplitudes and frequency's of excitation, the pendulum stays upright without falling over. Simulations provided below model this behavior and display the proper initial conditions that promote stabilization of the inverted mass. Animations are created to display the inverted pendulum stabilizing across variations in initial conditions.

Simulation Setup¶

Below can be seen two images that setup the problem we are looking at. The diagram of the model can be seen with the excited base moving up the Y axis and inverted mass connected with a rod free to move in the X-Y plane. The mass is only influenced by the force of gravity and restricted by the connection to the exicted base. The angle $\theta$ is defined as positive from the vertical Y-axis moving clockwise in the diagram.

diagram

The bond graph below shows the physical model of the system and how each part is interacting such that we can derive equations to solve for the motion of the inverted mass.

bond graph

Equations used for simulations¶

Velocity equations for X, Y, and $\theta$ directions

$$v_y = \dot{y} - l \dot{\theta} \sin(\theta) $$$$v_x = l \dot{\theta} \cos(\theta) $$$$ \dot{\theta} = \frac{p_x}{l\cos(\theta)\ m} $$

Momentum equations

$$p_x = m l \dot{\theta} \cos(\theta) $$$$p_y = m\left[\dot{y} - \frac{\sin(\theta)\ p_x}{\cos(\theta)\, m}\right] $$

Derivatives of momentum equations $$\dot{p_y} = m\left[\ddot{y} - \frac{\sin(\theta)\ \dot{p_x}}{\cos(\theta)\ m} - \frac{p_x\ \dot{\theta}}{m\ \cos^2{\theta}} \right] $$ $$\dot{p_x} = m g \sin(\theta)\cos(\theta) + m \ddot{y}\sin(\theta)\cos(\theta) - \frac{\sin(\theta)}{\cos(\theta)} \dot{\theta} p_x$$

Code Setup¶

Since we plan to do several simulations and visualize them, in this section we define several functions that allow us to simplify how this is done.

Define Derivative Function¶

To solve this non-linear system of ordinary differential equations, we need to start from an initial state, numerically integrate along the state derivatives to find the next state, and repeat as desired

For integration, we will use scipy.integrate's solve_ivp(), which takes in a range of time (t_span) through which to integrate, a function handle (func) that takes in the time and state and returns the derivatives, and the initial condition (initial) of the state.

Since we wish to solve this problem using several different parameters, we will define a function get_func() that returns a function handle that we can pass to solve_ivp(). We also return a function that includes additional state information.

In [1]:
from copy import deepcopy
import numpy as np

def get_func(params: dict[str, float]):
    # Make a copy of the input parameters to ensure the returned function 
    # does not change when the original parameters dictionary changes
    params = deepcopy(params)

    # Extract parameters from dict
    theta = params["theta"]
    ampl = params["amplitude"]
    freq = params["frequency"]*2*np.pi # Convert from Hz to rad/s
    leng = params["l_arm"]
    mass = params["m_mass"]
    g = params["g"]

    # Get initial condition (add position into the system so solve_ivp can integrate it for us)
    px = 0
    initial = [px, theta]

    # Create derivative calculator function to be used by solve_ivp()
    # state = [px, theta]
    def func(t, state):
        state = deepcopy(state)
        px = state[0]
        theta = state[1]

        # Sine and cosine of theta, for computational speed
        sin_theta = np.sin(theta)
        cos_theta = np.cos(theta)
        tan_theta = np.tan(theta)

        # Position of cart and its derivatives
        X = 0
        Y = ampl*np.sin(freq*t)
        d_Y = ampl*freq*np.cos(freq*t)
        dd_Y = -(freq**2)*Y

        # State derivatives
        d_theta = px/(leng*cos_theta*mass)
        d_px = mass*g*sin_theta*cos_theta + mass*dd_Y*sin_theta*cos_theta - tan_theta*d_theta*px

        # Other state variables
        x = X + leng*np.sin(theta) # x-position of pendulum mass
        y = Y + leng*np.cos(theta) # y-position of pendulum mass
        py = mass*(d_Y - tan_theta*px/mass) # y-momentum of pendulum mass
        d_py = mass*(dd_Y - tan_theta*d_px/mass - px*d_theta/(mass*cos_theta**2)) # first time-derivative of y-momentum of pendulum mass
        
        # Concatenate state derivatives
        d_state = [d_px, d_theta]

        # Create current state dict
        state = {}
        state["t"] = t
        state["X"] = X
        state["Y"] = Y
        state["x"] = x
        state["y"] = y
        state["px"] = px
        state["py"] = py
        state["theta"] = theta
        state["d_Y"] = d_Y
        state["d_px"] = d_px
        state["d_py"] = d_py
        state["d_theta"] = d_theta
        state["dd_Y"] = dd_Y
        state["tan_theta"] = tan_theta

        return d_state, state

    def func_wrap(t, state):
        '''
        Since solve_ivp() needs a function that only returns the state
        derivatives, we create a wrapper for func() that discards the rest of
        the state and returns only the state derivatives
        '''
        d_state, state = func(t, state)
        return d_state

    return func, func_wrap, initial
    

Propagate Solution Using solve_ivp()¶

Now that we have a function that calculates the state derivatives given the state, we need to propagate the solution throughout time. The function solve_problem() was designed to use the get_func() function we defined earlier, calculate the solution using scipy.integrate's solve_ivp(), and consolidate the interesting parts of the solution and the parameters given to define it.

In [2]:
from scipy.integrate import solve_ivp
import pandas as pd
import numpy as np


def solve_problem(
        params: dict[str, float],
        t_eval: np.ndarray,
        name: str,
        rtol: float = 1e-12,
        atol: float = 1e-12,
    ):

    func, func_wrap, initial = get_func(params)
    output = solve_ivp(
        func_wrap,
        (min(t_eval), max(t_eval)),
        initial,
        t_eval=t_eval,
        method="RK45",
        rtol=rtol,
        atol=atol,
    )

    ts = output.t
    ys = output.y

    states: list[dict[str, float]] = []
    for t, y in zip(ts, ys.T):
        d_state, state = func(t, y)
        states.append(state)

    df = pd.DataFrame(states)

    solution = {
        "params": params.copy(),
        "data": df,
        "name": name,
    }

    return solution

Define Plotting Functions¶

To visualize the solutions, we can create a static plot or an animation of the pendulum. Since we have several solutions we would like to plot together, we create functions to consolidate the code. plot_time() and plot_ani() each take in a list of dictionaries outputted as the solution from solve_problem().

In [3]:
from typing import Any
from matplotlib.lines import Line2D
import matplotlib.pyplot as plt
from matplotlib import animation
import numpy as np
import pandas as pd


def plot_time(solutions):
    fig = plt.figure()
    ax = fig.add_subplot()
    for solution in solutions:
        name = solution["name"]
        df = solution["data"]

        t_vals = df.get("t").to_numpy()
        theta_vals = df.get("theta").to_numpy() * 180 / np.pi

        ax.plot(t_vals, theta_vals, label=f"theta: {name}")

    plt.legend()
    plt.title("Theta vs. Time")
    plt.xlabel("time (s)")
    plt.ylabel("theta (deg)")
    return fig


def plot_ani(solutions: list[dict[str, Any]], interval: int = 10):
    fig = plt.figure()
    ax = fig.add_subplot()
    all_x = np.array([])
    all_y = np.array([])

    # Initialize arrays to store the animation plots and data to use to update them
    num_solutions = len(solutions)
    plots: list[dict[str, Line2D]] = [{} for n in range(num_solutions)]
    data: list[dict[str, np.ndarray]] = [{} for n in range(num_solutions)]

    for i, solution in enumerate(solutions):
        name: str = solution["name"]
        df: pd.DataFrame = solution["data"]

        data[i]["t_vals"] = df.get("t").to_numpy()
        data[i]["X_vals"] = df.get("X").to_numpy()
        data[i]["Y_vals"] = df.get("Y").to_numpy()
        data[i]["x_vals"] = df.get("x").to_numpy()
        data[i]["y_vals"] = df.get("y").to_numpy()

        (plots[i]["cart"],) = ax.plot([], [], "o", c="b", label=f"cart: {name}")
        (plots[i]["line"],) = ax.plot([], [], "-", c="k")
        (plots[i]["pend"],) = ax.plot([], [], "o", label=f"pendulum: {name}")

        num_frames = len(data[i]["t_vals"])

        all_x = np.concatenate((all_x, data[i]["X_vals"].copy(), data[i]["x_vals"].copy()))
        all_y = np.concatenate((all_y, data[i]["Y_vals"].copy(), data[i]["y_vals"].copy()))

    def update_points(n):
        returns = []
        for i in range(num_solutions):
            plots[i]["cart"].set_data(([data[i]["X_vals"][n]], [data[i]["Y_vals"][n]]))

            plots[i]["line"].set_data(
                (
                    [data[i]["X_vals"][n], data[i]["x_vals"][n]],
                    [data[i]["Y_vals"][n], data[i]["y_vals"][n]],
                )
            )

            plots[i]["pend"].set_data(([data[i]["x_vals"][n]], [data[i]["y_vals"][n]]))
            returns.extend([plots[i]["cart"], plots[i]["line"], plots[i]["pend"]])

        returns = tuple(returns)
        return (*returns,)

    ani = animation.FuncAnimation(
        fig, update_points, num_frames, interval=interval, blit=True, repeat=True
    )

    min_x, max_x = min(all_x), max(all_x)
    min_y, max_y = min(all_y), max(all_y)
    plt.xlim(min_x - 0.5, max_x + 0.5)
    plt.ylim(min_y - 0.5, max_y + 0.5)
    plt.legend()
    plt.title("Pendulum Animation")
    plt.xlabel("x (m)")
    plt.ylabel("y (m)")
    return ani

Simulation¶

Validation (Unexcited Pendulum)¶

Before running the simulation with the cart moving, we can run it without the moving cart to see if the pendulum behaves as expected, neither gaining nor losing energy.

In [4]:
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams["animation.html"] = "jshtml"

# -----------------------------
# Global Parameters
# -----------------------------
params: dict[str, float] = {}
params["l_arm"] = 1.0   # length of pendulum arm [m]
params["m_mass"] = 1.0  # mass of pendulum mass [kg]
params["g"] = 9.8       # gravitational force [m/s^2]

# Time
t_start: float = 0
t_end: float = 25
t_increment: float = 0.05 # visualization time step duration (sec)
# NOTE: playback speed is t_increment/(interval/1000)

t_span = (t_start, t_end)
t_eval = np.arange(min(t_span), max(t_span)+t_increment, t_increment)

params["amplitude"] = 0
params["frequency"] = 0
params["theta"] = np.deg2rad(60)

solution = solve_problem(params, t_eval, "validation")

fig = plot_time([solution])

As we can see from the validation case, the pendulum is behaving as expected. It does not appear to have any unknown forces acting on it, and it has a minimum and maximum angle of +/- 60°, which is where it was released.

10°, 20°, and 30° With Oscillating Cart¶

Here we simulate the inverted pendulum with the oscillating cart at the recommended amplitude and frequency.

In [5]:
# Time
t_start: float = 0
t_end: float = 5
t_increment: float = 0.005 # visualization time step duration (sec)
# NOTE: playback speed is t_increment/(interval/1000)

t_span = (t_start, t_end)
t_eval = np.arange(min(t_span), max(t_span)+t_increment, t_increment)


# Run simulations
solutions = []
params["amplitude"] = 0.1
params["frequency"] = 20
params["theta"] = np.deg2rad(10)
solutions.append(solve_problem(params, t_eval, "10°"))
params["amplitude"] = 0.1
params["frequency"] = 20
params["theta"] = np.deg2rad(20)
solutions.append(solve_problem(params, t_eval, "20°"))
params["amplitude"] = 0.1
params["frequency"] = 20
params["theta"] = np.deg2rad(30)
solutions.append(solve_problem(params, t_eval, "30°"))

fig = plot_time(solutions)

We see that the for the three initial angles, 10°, 20°, and 30° away from vertical, the inverted pendulum is stable with amplitude 0.1 meter and frequency 20 Hz.

Testing 40°¶

Although not necessary, it is interesting to see where the 0.1 meter, 20 Hz oscillation stops working. Adding in a 40° initial angle to the simulation shows that the oscillation no longer works at this initial angle. Here, we see the 40° pendulum begins to rotate continuously around the cart, as the cart maintains a higher level of energy in the system than it starts with.

In the case of 10°, 20°, and 30°, we see a sort of "destructive interference" between the oscillation of the pendulum and the oscillation of the cart. However, when we get to 40°, we see "constructive interference". This implies that there would be a better frequency/amplitude for a pendulum starting at this initial angle.

In [6]:
# Time
t_start: float = 0
t_end: float = 1.5
t_increment: float = 0.005 # visualization time step duration (sec)
# NOTE: playback speed is t_increment/(interval/1000)

t_span = (t_start, t_end)
t_eval = np.arange(min(t_span), max(t_span)+t_increment, t_increment)


# Run simulations
solutions = []
params["amplitude"] = 0.1
params["frequency"] = 20
params["theta"] = np.deg2rad(10)
solutions.append(solve_problem(params, t_eval, "10°"))
params["amplitude"] = 0.1
params["frequency"] = 20
params["theta"] = np.deg2rad(20)
solutions.append(solve_problem(params, t_eval, "20°"))
params["amplitude"] = 0.1
params["frequency"] = 20
params["theta"] = np.deg2rad(30)
solutions.append(solve_problem(params, t_eval, "30°"))
params["amplitude"] = 50
params["frequency"] = 0.1
params["theta"] = np.deg2rad(40)
solutions.append(solve_problem(params, t_eval, "40°"))

fig = plot_time(solutions)

Just for fun, we can use our animate function to show a time-based animation of all four pendulums. Before starting it, note the starting positions are prescribed, and then see the 40° pendulum make full rotations about the cart.

In [7]:
from IPython.display import HTML

# Plot solutions
ani = plot_ani(solutions, interval = 15)
plt.axis('equal')

HTML(ani.to_jshtml())
Out[7]:

Frequency/Amplitude Space Exploration for 40°¶

While the 0.1 meter, 20 Hz oscillation did not make the 40° cart stable in the inverted position, there might be a combination of amplitude/frequency that it works at. Here, we sweep a frequency and amplitude range to see which if any combination might work.

In [ ]:
from concurrent.futures import ThreadPoolExecutor
from concurrent.futures import as_completed
from itertools import product
from tqdm import tqdm

# Time
t_start: float = 0
t_end: float = 1.5
t_increment: float = 0.020 # visualization time step duration (sec)
# NOTE: playback speed is t_increment/(interval/1000)

t_span = (t_start, t_end)
t_eval = np.arange(min(t_span), max(t_span)+t_increment, t_increment)

num_frequencies = 20
num_amplitudes = 20
frequencies = np.linspace(0.01, 100, num_frequencies)
amplitudes = np.linspace(0.01, 100, num_amplitudes)

F, A = np.meshgrid(frequencies, amplitudes)
T = np.zeros_like(A)

pbar = tqdm(total=num_frequencies*num_amplitudes)

def process(params: dict, i, j, freq, ampl):
    params.copy()
    params["frequency"] = freq
    params["amplitude"] = ampl
    solution = solve_problem(params, t_eval, "sweep")
    result = np.max(np.abs(solution["data"]["theta"]))
    pbar.update()
    return result

params["theta"] = np.deg2rad(40)
results = {}
with ThreadPoolExecutor(max_workers=2) as executor:
    for j, i in product(range(num_frequencies), range(num_amplitudes)):
        freq = F[i, j]
        ampl = A[i, j]

        results[executor.submit(process, params, i, j, freq, ampl)] = [i, j]

    for future in as_completed(results):
        i, j = results[future]
        data = future.result()

T_min, T_max = T.min(), T.max()

fig, ax = plt.subplots()

c = ax.pcolormesh(F, A, T, cmap='RdBu', vmin=T_min, vmax=T_max)
# set the limits of the plot to the limits of the data
ax.axis([F.min(), F.max(), A.min(), A.max()])
ax.set_title('Theta Maximum Magnitude')
ax.set_xlabel("Frequency (Hz)")
ax.set_ylabel("Amplitude (m)")
fig.colorbar(c, ax=ax)

plt.show()
  0%|          | 0/400 [00:00<?, ?it/s]  9%|▉         | 37/400 [02:15<57:36,  9.52s/it]  

As we can see from the heatmap, there are several possible combinations of amplitude and frequency that satisfy the inverted pendulum stability requirement. Dark red values are combinations that keep the maximum angle of the pendulum less than 90°.

Results and Discussion¶

The results from the simulations above describe the behavior that we would be expecting from the model created. We can see that for set initial angles of $\theta$ = 10, 20, and 30 degrees we were able to recover the inverted pendulum to stay upright. Starting at 90 was too far from vertical that no amplitude or excitation was able to recover it. Varying numbers of amplitude and frequency of the excitation.

Conclusions¶

The inverted pendulum simulation is a great example of behavior that seems non-physical but can be a reality given the right set-up conditions to make it possible. It is also a great test case for bond graph modelling that enables simpler descriptions of the system to be simulated. As can be seen, the proper modelling and simulation of the system allows use to view the behavior of the pendulum for different starting conditions and how it has the possibility to stay upright.